
#reading in fits data file
from astropy.coordinates import SkyCoord
from astropy.io import fits 
import matplotlib.pyplot as plt
from astropy.table import Table
import numpy as np


#reading in fits files 

#wavelength solutions 
wavedata= fits.getdata('/Users/sandratredinnick/Desktop/ASTR5550/z2r50205t.c0f.fits')
print(wavedata)
headerwave= fits.getheader('/Users/sandratredinnick/Desktop/ASTR5550/z2r50205t.c0f.fits')
print(wavedata.shape)

#calibrated fluxes
calfluxdata= fits.getdata('/Users/sandratredinnick/Desktop/ASTR5550/z2r50205t.c1f.fits')
print(calfluxdata)

#uncertainities
uncertdata= fits.getdata('/Users/sandratredinnick/Desktop/ASTR5550/z2r50205t.c2f.fits')
print(uncertdata)

[[2792.38625404 2792.39387204 2792.40148979 ... 2807.30055332
  2807.3078941  2807.31523656]
 [2792.38625404 2792.39387204 2792.40148979 ... 2807.30055332
  2807.3078941  2807.31523656]
 [2792.38625404 2792.39387204 2792.40148979 ... 2807.30055332
  2807.3078941  2807.31523656]
 ...
 [2791.96152474 2791.96914634 2791.97676771 ... 2806.88311864
  2806.89046314 2806.89780933]
 [2791.96152474 2791.96914634 2791.97676771 ... 2806.88311864
  2806.89046314 2806.89780933]
 [2791.96152474 2791.96914634 2791.97676771 ... 2806.88311864
  2806.89046314 2806.89780933]]
(16, 2000)
[[2.2250016e-13 3.0200262e-13 6.3603928e-14 ... 1.7262567e-13
  5.5018222e-14 2.1968450e-13]
 [1.4303592e-13 1.4305775e-13 1.4307957e-13 ... 3.3727434e-13
  2.6671747e-13 2.4320300e-13]
 [3.7150437e-13 2.9208875e-13 5.3690247e-14 ... 1.8643812e-13
  1.3939940e-13 2.8054525e-13]
 ...
 [1.3972841e-13 2.1924210e-13 1.3977200e-13 ... 2.6314199e-13
  9.8672908e-14 1.9266676e-13]
 [1.3913449e-13 2.1864817e-13 2.1868071e-13 ... 1.9207281e-13
  1.6858249e-13 1.9208588e-13]
 [3.7853298e-13 2.9909487e-13 1.4013291e-13 ... 9.9022244e-14
  1.6951636e-13 1.9301975e-13]]
[[1.37627336e-13 1.58939851e-13 7.94806658e-14 ... 6.65300876e-14
  4.07424779e-14 7.43876261e-14]
 [1.12372263e-13 1.12387462e-13 1.12402648e-13 ... 9.11000740e-14
  8.14849558e-14 7.80183955e-14]
 [1.77676152e-13 1.58939851e-13 7.94806658e-14 ... 7.05658134e-14
  6.22351698e-14 8.48149404e-14]
 ...
 [1.12403434e-13 1.37683986e-13 1.12433649e-13 ... 8.13928935e-14
  5.25405186e-14 7.04926704e-14]
 [1.12403434e-13 1.37683986e-13 1.37702539e-13 ... 7.04883200e-14
  6.64590724e-14 7.04926704e-14]
 [1.77725429e-13 1.58983774e-13 1.12433649e-13 ... 5.25388957e-14
  6.64590724e-14 7.04926704e-14]]


#plot data first to see what im working with

#plot first row of wavelengths vs flux 
##16 rows so it seems like those are my 16 sub exporesure 
plt.figure(figsize=(20,10), dpi=400)
plt.plot(wavedata[0],calfluxdata[0])
plt.xlabel("Wavelength")
plt.ylabel("Flux")

Text(0, 0.5, 'Flux')


plt.figure(figsize=(20,10), dpi=400)
plt.plot(wavedata[1],calfluxdata[1])
plt.xlabel("Wavelength")
plt.ylabel("Flux")

Text(0, 0.5, 'Flux')


plt.figure(figsize=(20,10), dpi=400)
plt.plot(wavedata[2],calfluxdata[2])
plt.xlabel("Wavelength")
plt.ylabel("Flux")

Text(0, 0.5, 'Flux')


plt.figure(figsize=(20,10), dpi=400)
plt.plot(wavedata[3],calfluxdata[3])
plt.xlabel("Wavelength")
plt.ylabel("Flux")

Text(0, 0.5, 'Flux')


#finding the 4 groups of the same wavelenth (so i can co-add spectra)
## co-add: mean of the flux arrays 0-3,4-7,8-11,12-15
print(wavedata[0]==wavedata[3])
print(wavedata[4]==wavedata[7])
print(wavedata[8]==wavedata[11])
print(wavedata[12]==wavedata[15])

[ True  True  True ...  True  True  True]
[ True  True  True ...  True  True  True]
[ True  True  True ...  True  True  True]
[ True  True  True ...  True  True  True]


#wavedata[0]==wavedata[2]


#make array of 4 sets to co add spectra
array1= np.array(calfluxdata[0:4])
array2= np.array(calfluxdata[4:8])
array3=np.array(calfluxdata[8:12])
array4= np.array(calfluxdata[12:16])

#co adding spectra 
setone= np.mean(array1,axis=0)
settwo= np.mean(array2,axis=0)
setthree= np.mean(array3,axis=0)
setfour= np.mean(array4,axis=0)


#n= index 
#f = flux array one 
#g = flux array two 
#sum over all values in m
#shift by m


####GHRS.pdf document said:The best way to align the subexposures
#is by cross-correlating them 
#against the first exposure obtained

#create cross correlation function because cant use np or scipy correlate funct.
N=100

def crosscorr(spectra_one,spectra_two):
    crosscorr=np.zeros(N)
    for i in range (N):
        x= np.roll(spectra_two,-i) #np.roll inverts the second array flux 
        y=x*spectra_one #multiply each element above by element of other flux array 
        crosscorr[i]=sum(y) #sum over the multiplication 
    return crosscorr


#cross correlate the fluxes against the first exposure obtained 
ccfone= crosscorr(setone,settwo)
ccftwo=crosscorr(setone,setthree)
ccfthree= crosscorr(setone,setfour)


#shift by the number of wavelength bins 
#required to achieve maximum cross-correlation

#find the shift, number of wavelength binds required to achieve max correlation

shiftone = np.argmax(ccfone)
shifttwo = np.argmax(ccftwo)
shiftthree = np.argmax(ccfthree)
print(shiftone)
print(shifttwo)
print(shiftthree)

17
33
52


# allign/shifting subexposures 
allignone= np.roll(settwo,-shiftone)
alligntwo= np.roll(setthree,-shifttwo)
allignthree= np.roll(setfour,-shiftthree)


#co adding all alligned subexposures 
arraycoadd=np.array([setone,allignone,alligntwo,allignthree])
coadd2= np.mean(arraycoadd,axis=0) #this is my master flux array


#calculating uncertainities in co adding 
#need to use propogation of errors --> boy am I am getting good at prop of errors! 
unccoaddone = ((uncertdata[0])**2+(uncertdata[1])**2+(uncertdata[2])**2+(uncertdata[3])**2)/16
#print(unccoaddone)
unccoaddtwo = (uncertdata[4]**2+uncertdata[5]**2+uncertdata[6]**2+uncertdata[7]**2)/16
#print(unccoaddtwo)
unccoaddthree = (uncertdata[8]**2+uncertdata[9]**2+uncertdata[10]**2+uncertdata[11]**2)/16
#print(unccoaddthree)
unccoaddfour = (uncertdata[12]**2+uncertdata[13]**2+uncertdata[14]**2+uncertdata[15]**2)/16
#print(unccoaddfour)


# shifting 
unccoaddshiftone = np.roll(unccoaddtwo, -shiftone)
unccoaddshifttwo = np.roll(unccoaddthree, -shifttwo)
unccoaddshiftthree = np.roll(unccoaddfour, -shiftthree)


#finalized calculated uncertainity values 
sigma2_unccoaddshift= (unccoaddone+ unccoaddshiftone+ unccoaddshifttwo+unccoaddshiftthree)/16
sigma_unccoaddshift = np.sqrt(sigma2_unccoaddshift)


#plot co added spectrum with uncertainities 
plt.figure(dpi = 200, figsize = (8, 6))
plt.plot(wavedata[0],coadd2,label="Spectrum")
plt.errorbar(wavedata[0],sigma_unccoaddshift,label="Uncertainity")
plt.fill_between(wavedata[0],coadd2-sigma_unccoaddshift,coadd2+sigma_unccoaddshift,color='black',label="Uncertainity")
plt.title("Co-added Spectrum With Associated Uncertainity Values")
plt.xlabel("Wavelength ($\AA$)")
plt.ylabel("Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)")
plt.legend()

<matplotlib.legend.Legend at 0x7fe278a01310>


#find s/n 
#s = master flux array 
#n= errors

#recall: coadd2 is the alligned co added fluxes 
master_flux = coadd2

signoise= master_flux/sigma_unccoaddshift
print(signoise)
print(np.mean(signoise))

[ 9.59806   9.402988  9.870239 ... 12.560999 12.236406 13.652048]
12.339193


#plot s/n 
plt.figure(dpi = 200, figsize = (8, 6))
plt.plot(wavedata[0],signoise)
plt.xlabel("Wavelength ($\AA$)")
plt.ylabel("Signal to Noise Ratio (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)")
plt.title("Spectrum of Signal to Noise Ratio as a Function of Wavelength")
#plt.axvline(2796.3, color = 'black', label='Mg II absorption line')
#plt.axvline(2803.48, color = 'black', label='Mg II absorption line')

Text(0.5, 1.0, 'Spectrum of Signal to Noise Ratio as a Function of Wavelength')


#need to find indicies of absoption lines! 
#note the absoption lines are a range, not just 
#a single horizontal line 
np.argmin(signoise),1000+ np.argmin(signoise[1000:])

(515, 1479)


# find the average continuum flux value between 
#the two Mg II lines, and divided by the standard deviation in that region

#find the signal to nosie array between the abs lines 
#so did 515+5, 1479-5 to get inbetween array 
signoisebetween = signoise[520:1474]


#mean contiuum value beetween the 2 abs lines 
meansignoisebetween = np.mean(signoisebetween)
print(meansignoisebetween)

#std vaue of between the 2 abs line contiuum 
stdsignoisebetween= np.std(signoisebetween)
print(stdsignoisebetween)

12.324845
0.61785144


#find s/n between the two absorption lines for flux 
#not using 515:1480 bc the abs line is a range of values not 
#just a single point
fluxbetween= master_flux[520:1474]

#mean of the the between contiuum 
meanfluxbetween= np.mean(fluxbetween)
print(meanfluxbetween)


# find the standard devition of the between continuum 
stdfluxbetween= np.std(master_flux[520:1474])
print(stdfluxbetween)


#find the s/n for the contiuum between the 2 abs lines 
print(meanfluxbetween/stdfluxbetween)

2.326488e-13
1.9952174e-14
11.660323


percent = (abs((meanfluxbetween-3*stdfluxbetween-meanfluxbetween)/meanfluxbetween))*100
print(percent)

25.728274398227423


#find s/n contintuum (array) on either side of the absoption lines 

#goal: find s/n on left and right sides of absoption lines
#so need mean and standard deviation of left and right sides! (s/n)

#recall the absoption line is a range of values not a single point
#first--> left side 
#find array that is the left side of asborption line 
spectrumleft= coadd2[:510]
#finding mean of spectra to the left side of absoption line
spectrumleftmean = np.mean(spectrumleft)
#findind std of spectra to the left side of absoption line 
spectrumleftstd = np.std(spectrumleft)
#print(spectrumleftmean)
#print(spectrumleftstd)

#s/n for left side of asbotion lines (0 to 514)
spectrumleftSN = spectrumleftmean/spectrumleftstd
print(spectrumleftSN)


#second --> right side 
#find array that is the right side of asborption line 
spectrumright= coadd2[1480+5:]
#finding mean of spectra to the right side of absoption line
spectrumrightmean = np.mean(spectrumright)
#findind std of spectra to the right side of absoption line 
spectrumrightstd = np.std(spectrumright)
#print(spectrumrightmean)
#print(spectrumrightstd)

#s/n for right side of asbotion lines (1480 to end)
spectrumrightSN = spectrumrightmean/spectrumrightstd
print(spectrumrightSN)

12.267691
12.0922575


#find the S/N in two absorption line 

#first find the S/N in the left absoption line (line near 2796 ̊A)
#signoise = signal to noise of whole continum so just
#get s/n of the absoption line 
leftabslineSN = signoise[515]
print(leftabslineSN)

#second find the s/n for the right asboption line (line near 2803  ̊A)
rightabslineSN = signoise[1479]
print(rightabslineSN)

5.1924796
7.073334


#fractional uncertainity stuff

#goal: assuming the noise is dominated by Poisson noise from photon counts
##need to compute fractional uncertainities in the continuum and the line

#use def of fractional uncertanity, and to look at the relationship between them
#take the ratio --> fractional uncertainities in continuum/fractional unct in line
#note that uncertainities for poisson stats are defined by sqrt(mu)! 


#deciding to the continuum inbetween the absorption lines 
#calculating the ratio for the line, chooshing the left absorption line 
ratiocontline =  np.sqrt(coadd2[515])/np.sqrt(np.mean(coadd2[520:1469]))
print(ratiocontline)

0.47227946


#find s/n ratio for the conintuum inbetween the 2 abs lines 
SNratio = (np.std(coadd2[520:1469])*coadd2[515])/(np.mean(coadd2[520:1469])*sigma_unccoaddshift[515])
print(SNratio)

#now can compare!

0.44459233


#first step: mask out the absorption lines 


spectrumleftnew= np.where(wavedata[0]<wavedata[0][500])[0]

spectrumbetweennew= np.where((wavedata[0]>wavedata[0][530]) & (wavedata[0]<wavedata[0][1454]))[0]

spectrumrightnew = np.where(wavedata[0] > wavedata[0][1494])[0]

indicies = [spectrumleftnew,spectrumbetweennew,spectrumrightnew]
indicies2 = np.concatenate(indicies)

meanspectrummask= coadd2[indicies2]
wavesmask= wavedata[0][(indicies2)]


#plot the spectrum or continuum of the two absoprtion 
#lines masked out 
#note that the line inbetween is Python artifically connecting 
#the data points but the line doesnt have meaning 
plt.plot(wavesmask, meanspectrummask, color='black')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)')
plt.title('Spectrum With Asborption Lines Masked Out')
plt.show()


#doing a polynomial fit 
#doing a degree one polynomial fit 

#array for params and fit 
parameters,fit = np.polyfit(wavesmask,meanspectrummask,deg=1,cov=True)
#creating an array for the fit to make a smooth curve 
x = np.linspace(wavesmask[0],wavesmask[-1],1000)
#creating the line fit
##parameters[0]=slope 
y = parameters[0]*x + parameters[1]


#plot masked spectrum with fit overplotted 
plt.figure(dpi = 200, figsize = (8, 6))
plt.plot(wavesmask, meanspectrummask, color='black',label='Specturm')
plt.plot(x,y,label='Degree One Polynomial Fit')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)')
plt.legend()
#plt.axhline(np.mean(meanspectrummask), c='green',label='Mean of spectrum')
plt.title('Spectrum With Asborption Lines Masked Out & Polynomial Fit')
plt.show()


#can fix code to be different, take sqrt inside 
#uncertainities in fit paramters 

sig2slopeone = fit[0][0]
sig2inttwo = fit[1][1]
#uncertainity in slope paramater 
sigslopeone = np.sqrt(sig2slopeone)
#uncertainity in intercept paramater 
siginttwo = np.sqrt(sig2inttwo)
print(sigslopeone)
print(siginttwo)

#print paramter estimates
#slope
print(parameters[0])
#intercept
print(parameters[1])

1.0115119626338277e-16
2.832114711126018e-13
2.8235668811163627e-16
-5.583116058824341e-13


#creating best fit line for normalized specturm, abs lines
#masked out
bestfit = parameters[0]*wavesmask+parameters[1]

#normalizing the spectrum with abs lines masked out 
normspec= meanspectrummask/bestfit
normspec

#sanity checking to make sure it is normalized. 
#if it was normalized would get mean =1
meannorm = np.mean(normspec)
print(meannorm)

0.999999889775847


#plotting the normalized spectrum 
plt.figure(dpi = 200, figsize = (8, 6))
plt.plot(wavesmask, normspec, color='black')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)')
plt.title('Normalized Spectrum With Asborption Lines Masked Out')
#plt.axhline(meannorm, c='green',label='Mean of spectrum')
plt.show()


#goal: do gaussian fit to both absorption lines 

#step one: first need normalized spectrum with 
#the asbroption lines in 
###so I must renormalize with a best fit line that considers
#the whole spectrum instead of the specturm with the absroption 
#lines removed 

#same procedure as 2.4
#getting best fit line, using whole spectrum with 
#absorption lines in 
bestfit2 = parameters[0]*wavedata[0]+parameters[1]

#normalized spectrum array with abs lines in 
normspec2= coadd2/bestfit2

#sanity check--> make sure its actually normalized by checking mean
#mean should be close to 1!
meannorm2 = np.mean(normspec2)


#plot Normalized Spectrum/continuum with Abs Lines in spectrum 
plt.plot(wavedata[0], normspec2, color='black')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)')
plt.title('Normalized Spectrum With Absoption Lines')
#plt.axhline(meannorm, c='green',label='Mean of spectrum')
plt.show()


#now I need to find two different arrays where the two different arrays 
#are the first absoption line and second abs
#####similar process to 2.3, but flipped, want only the absoption lines 
#note I am using the same code as 2.3 to get these two arrays (the code i used
#from 2.3 was Conis and I am applying it again here, see 2.3 for my explanation 
#of exactly what the code is doing)


#first absorption line array
lineone= np.where((wavedata[0]>wavedata[0][505]) & (wavedata[0]<wavedata[0][525]))[0]
print(lineone)

#first absoption line flux array and wavelength array 
##note that I multiplied the flux array by -1 because I need to invert the flux array 
#to properly fit a gaussian 
#fluxarray
fluxlineone= -normspec2[lineone]+1 #shift up by one to get on right scale 
#wavelength array 
waveslineone= wavedata[0][lineone]

[506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523
 524]


#plot the absoption line continuum for the first line (line to the left of the spectrum)

#plt.scatter(waveslineone, fluxlineone, color='black')
#plt.xlabel('Wavelength')
#plt.ylabel('Flux')
#plt.title('Spectrum of Absorption Line One (Near 2796 $\AA$)')
#plt.show()


from scipy.optimize import curve_fit


#gaussian fit for the absoption lines 
## I found this code from stack overflow! (see link directly above)
#parameters of fit are: central wavelength (mean),amplitude, and width (standard deviation)


#first define the gaussian function 
def Gauss(x, amp, mean, sigma):
    return amp * np.exp(-(x - mean)**2 / (2 * sigma**2))

#x= wavelength array of the first abs line 
#amp= amplitude


#need starting values for the gaussian fit of abs line one! 
#fit by eye 
#worked together with coni to set the starting values 
am= np.max(fluxlineone)-np.min(fluxlineone)
mean= waveslineone[np.argmin(fluxlineone)]
std=.09


#doing the fitting for absorption line 1
param,unc = curve_fit(Gauss, waveslineone, fluxlineone, p0=[np.max(fluxlineone)-np.min(fluxlineone),waveslineone[np.argmin(fluxlineone)],.09])
#want to create an array of wavelengths to make the gaussian fit smooth for abs line 1
x1= np.linspace(min(waveslineone),max(waveslineone),1000)


#plotting the asboption line one continuum/spectrum with fit 
plt.plot(x1,-Gauss(x1,param[0],param[1],param[2])+1,c='red') #invert the gaussian fit
#function wont fit a inverse gaussian! sad! but working around it. 
plt.scatter(waveslineone, -fluxlineone+1, color='black') #invert the fluxes 
plt.title('Normalized Spectrum of Absorption Line One (Near 2796 $\AA$) With Gaussian Fit')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)')
plt.show()


#best fit values and uncertanities for fit parameters for abs line one 

#best fit values for amplitude, mean, standard deviation respectively
print(param)

#uncertainties associated with amplitude, mean, standard deviation respectively 
#the diagnoials are the uncertanties associated with each paramater 
print(unc)

[7.83477723e-01 2.79629393e+03 2.80868715e-02]
[[ 6.62321913e-04 -1.59543516e-07 -1.62450183e-05]
 [-1.59543516e-07  1.13232367e-06  1.49759747e-08]
 [-1.62450183e-05  1.49759747e-08  1.17286212e-06]]


#second absorption line array
linetwo= np.where((wavedata[0]>wavedata[0][1469]) & (wavedata[0]<wavedata[0][1489]))[0]
print(linetwo)

#first absoption line flux array and wavelength array 
##note that I multiplied the flux array by -1 because I need to invert the flux array 
#to properly fit a gaussian 
fluxlinetwo= -normspec2[linetwo]+1
waveslinetwo= wavedata[0][linetwo]

[1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483
 1484 1485 1486 1487 1488]


#plot the absoption line continuum for the second line (line to the left of the spectrum)
#plt.scatter(waveslinetwo, fluxlinetwo, color='black')
#plt.xlabel('Wavelength')
#plt.ylabel('Flux')
#plt.title('Spectrum of Absorption Line Two (Near 2803 $\AA$)')
#plt.show()


#need starting values for the gaussian fit for abs line 2 ! 
#fit by eye
#worked together with coni to set the starting values 
am3= np.max(fluxlinetwo)-np.min(fluxlinetwo)
mean3= waveslinetwo[np.argmin(fluxlinetwo)]
std3=.02


#doing the fitting for absorption line 2
param2,unc2 = curve_fit(Gauss, waveslinetwo, fluxlinetwo, p0=[np.max(fluxlinetwo)-np.min(fluxlinetwo),mean3,std3])

#want to create an array of wavelengths to make the gaussian fit smooth for abs line 1
x2= np.linspace(min(waveslinetwo),max(waveslinetwo),1000)


#plotting the asboption line two continuum/spectrum with fit 
plt.plot(x2,-Gauss(x2,param2[0],param2[1],param2[2])+1,c='red') #invert the gaussian fit
#function wont fit a inverse gaussian! sad! but working around it. 
plt.scatter(waveslinetwo, -fluxlinetwo+1, color='black') #invert the fluxes 
plt.title('Normalized Spectrum of Absorption Line Two (Near 2803 $\AA$) With Gaussian Fit')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)')
plt.show()


#goal: calculate the central velocity for each line 
#will need to use propogation of errors 
#best fit line one
###amp,mean,std=[7.83477723e-01 2.79629393e+03 2.80868715e-02]
##uncetainiies: 6.62321913e-04, 1.13232367e-06, 1.17286212e-06]

#best fit line two
##amp,mean,std= [6.19312544e-01, 2.80347993e+03, 2.49001433e-02]
##uncertainities: 1.60749018e-03, 3.45716898e-06, 3.48693944e-06

#observed = parameter values found in 3.1
#central velcoties for abs line one 
lambdarestone=  2796.352 #found in 2.1 
light= 2.99e+8 #speed of light 
centvelone= light*(param[1]-lambdarestone)/lambdarestone # in m/s
#errors for velecoties for abs lines one
##use prop of errors
uncvelone = (light/lambdarestone)*np.sqrt(unc[1][1])


#observed = parameter values found in 3.1
#central velcoties for abs line two 
lambdaresttwo= 2803.530 #found in 2.1 
centveltwo= light*(param2[1]-lambdaresttwo)/lambdaresttwo 
#errors for velecoties for abs lines 2 
##use prop of errors
uncveltwo = (light/lambdaresttwo)*np.sqrt(unc2[1][1])


print(centvelone)
print(uncvelone)
print(centveltwo)
print(uncveltwo)

-6209.12674621733
113.77966292471372
-5340.074961801732
198.30163008107883


#Test the consistency by using a differenced χ2 consistency test, 
#taking into account the uncertainties in the laboratory wavelengths of these lines

#to do the chisquared test we need the variance of g
##var_g= var_velocityabsone +var_velocityabstwo
#use propogation of errors 

#variance of g
sig2_g= ((light/lambdarestone)**2)*unc[1][1] + ((light*param[1]/lambdarestone**2)**2)*(0.001)**2 +((light/lambdaresttwo)**2)*unc2[1][1] + ((light*param2[1]/lambdaresttwo**2)**2)*(0.001)**2


#using notes from lecture 18 on canvas 
#finding the t statistic!
#assume mu is zero because doing difference in means test 
tstat = (centvelone-centveltwo)**2/sig2_g
print(tstat)

from scipy.stats import chi2
#fidning PTE for t Test
#DOF is N-1, N= 2 because have two velocity points 
DOF = 1 
PTE = 1 - (chi2.cdf(tstat,DOF))
print(PTE)

10.059830398142388
0.001515371760480999

Data Project 2: UV Spectral Line Analysis¶

Question One: Constructing a Spectrum for Analysis¶

1.1¶

1.2¶

1.3¶

1.4¶

Question Two: Normalizing by the Continuum¶

2.1¶

2.2¶

2.3¶

2.4¶

Question Three: Line Analysis¶

3.1¶

3.2¶